function moments = analyze_sensitivity(prefix_pf, prefix_skill, ...
                                       prefix, knitro_path, relative)

    % if boolean `relative` is true, relative changes to
    % parameters are used, otherwise absolute changes are used
    
    addpath('../')
    addpath('../Densities')
    addpath('../UtilityFunctions')
    addpath('../ProductionFunctions')
    addpath('../Economies/')
    addpath(knitro_path)

    dir_in = '../../../data/generated/matlab/Calibration/'
    dir_out = '../../../data/generated/matlab/simulation/'

    w0_range = 1e-4;
    wf_range = 1e4;

    % load calibrated parameters for skill distribution
    par_cal = tools.csv2struct([dir_in,prefix_skill,'_abilityDistThreeFixC_sol.csv']);
    par_mom = tools.csv2struct([dir_in,prefix_skill,'_abilityDistThreeFixC_mom.csv']);
    inequ_par = tools.csv2struct([dir_in,'inequ_par.csv']);
    constants = tools.csv2struct('../../../data/raw/constants.csv');

    DensSkill = abilityDistFixed(par_cal);

    DensPart = ParticipationNormal(par_cal);

    % utility function object
    par_util.epsilon = constants.epsilon;
    Util = UtilityQuasiLinear(par_util);

    % combine in SkillAndParticipation object
    par_dist.dens_skill = DensSkill;
    par_dist.dens_part = DensPart;
    par_dist.util = Util;

    par_dist.num_nodes_phi = 50;
    par_dist.scale_inc = par_cal.scale_inc; 

    par_dist.hsv_tax = true;
    Dist = SkillAndParticipation(par_dist);

    % Production function
    par_pf =  tools.csv2struct([dir_in,prefix_pf,'_KrusellRobotsEach_sol.csv'])
    mom_pf =  tools.csv2struct([dir_in,prefix_pf,'_KrusellRobotsEach_mom.csv'])

    par_pf.parse_inputs = true;
    par_pf.KrusellRobotsEach_use_properties = true;
    ProdFun = KrusellRobotsEach(par_pf);

    % Economy object
    par_econ.pf = ProdFun;
    par_econ.util = Util;

    par_econ.w_lb = 1e-6;
    par_econ.w_ub = 1e4;

    % robot prices and taxes
    inp.tau_B = par_pf.tau_B;
    inp.tau_E = par_pf.tau_E;
    inp.tau_S = par_pf.tau_S;

    par_econ.q_B = mom_pf.q_B;
    par_econ.q_E = mom_pf.q_E;
    par_econ.q_S = mom_pf.q_S;

    q_B_0 = par_econ.q_B;
    q_B_1 = mom_pf.Y_B_M_1/(1+inp.tau_B);

    par_econ.dist = Dist;

    par_econ.solver = 'knitro_eqsolve'

    par_econ.phi_lb = constants.phi_lb;
    par_econ.r = inequ_par.inequ_par;
    par_econ.par_dist = par_cal;

    par_econ.varnames = {'L_M','L_R','L_C', ...
                        'K_B_M','K_B_R','K_B_C', 'K_E','K_S',...
                        'tau_B','tau_E','tau_S','transfer','tau_hsv','lambda_hsv','compensation'};

    econ = KrusellWithRobotsEachEconomy(par_econ);
    econ.set_num_nodes(100,50,20,50);

    inp.L_M = par_pf.L_M_0;
    inp.L_R = par_pf.L_R_0;
    inp.L_C = par_pf.L_C_0;
    inp.K_B_M = par_pf.K_B_M_0;
    inp.K_B_R = par_pf.K_B_R_0;
    inp.K_B_C = par_pf.K_B_C_0;
    inp.K_E = par_pf.K_E_0;
    inp.K_S = par_pf.K_S_0;
    inp.transfer = par_cal.transfer;
    inp.tau_hsv = par_cal.tau_hsv;
    inp.lambda_hsv = par_cal.lambda_hsv;
    inp.compensation = 0;

    lb.L_M = 1e-16;
    lb.L_R = 1e-16;
    lb.L_C = 1e-16;
    lb.K_B_M = 1e-16;
    lb.K_B_R = 1e-16;
    lb.K_B_C = 1e-16;
    lb.K_E = 1e-16;
    lb.K_S = 1e-16;
    lb.tau_E = inp.tau_E;
    lb.tau_S = inp.tau_S;
    lb.transfer = inp.transfer;
    lb.tau_hsv = inp.tau_hsv;
    lb.lambda_hsv = inp.lambda_hsv;
    lb.compensation = 0;


    ub.L_M = 1e16;
    ub.L_R = 1e16;
    ub.L_C = 1e16;
    ub.K_B_M = 1e16;
    ub.K_B_R = 1e16;
    ub.K_B_C = 1e16;
    ub.K_E = 1e16;
    ub.K_S = 1e16;
    ub.tau_E = inp.tau_E;
    ub.tau_S = inp.tau_S;
    ub.transfer = inp.transfer;
    ub.tau_hsv = inp.tau_hsv;
    ub.lambda_hsv = inp.lambda_hsv;
    ub.compensation = 0;


    %% vesion with fixed capital (as sanity check for the wage-decile elasticities)
    inp.tau_B = par_pf.tau_B;
    lb.tau_B = inp.tau_B;
    ub.tau_B = inp.tau_B;

    droprange = sort(linspace(0,1-q_B_1/q_B_0,10));
    econ.fix_other_capital = true;
    
    % compute baseline moments (corresponds to current calibration)
    moments_baseline = compute_moments(econ, inp, lb, ub, q_B_0, q_B_1, droprange,w0_range,wf_range,prefix,'simulation_fix_capital',dir_out,par_cal, par_pf);
    moments= table();
    pars_to_vary = {'alpha',...
               'rho',...
               'sigma',...
               'kappa_M',...
               'kappa_R',...
               'kappa_C',...
               'nu',...
               'lambda',...
               'mu',...
               'eta_M',...
               'eta_R',...
               'eta_C',...
               'upsilon'};
           
    for par = pars_to_vary

        moments_tmp = adjust_parameter(par{1}, ...
            moments_baseline, ...
            econ, ...
            inp, ...
            lb, ...
            ub, ...
            q_B_0, ...
            q_B_1, ...
            droprange, ...
            w0_range, ...
            wf_range, ...
            prefix, ...
            'simulation_fix_capital', ...
            dir_out, ...
            par_cal, ...
            par_pf, ...
            relative);
        moments = [moments; moments_tmp];
 
    end
    
    % use translated parameter names (corresponding to the ones
    % used in the paper)
    pars_to_vary_translated = {'alpha',...
               'rho',...
               'sigma',...
               'kappa_M',...
               'kappa_R',...
               'kappa_C',...
               'kappa_MR',...
               'gamma_E',...
               'gamma_C',...
               'gamma_B_M',...
               'gamma_B_R',...
               'gamma_B_C',...
               'gamma_M'};
           
    
    % add column with name of varied parameter
    moments.par_changed = string(pars_to_vary_translated');
    
    % write to disk
    moments_struct = table2struct(moments);
    tools.struct2csv(moments_struct,[dir_out,prefix,'_moments_sensitivity.csv'])
end

function moments = adjust_parameter(par_to_adjust, ...
    moments_baseline, ...
    econ, ...
    inp, ...
    lb, ...
    ub, ...
    q_B_0, ...
    q_B_1, ...
    droprange, ...
    w0_range, ...
    wf_range, ...
    prefix, ...
    suffix, ...
    dir_out, ...
    par_cal, ...
    par_pf, ...
    relative)

    par_pf_changed = par_pf; % reset to original parameters
    val_to_adjust = par_pf.(par_to_adjust);
    
    if relative
        % use relative adjustment, leading to elasticities
        adj_factor = 0.99;

        if ismember(par_to_adjust, {'mu', 'eta_M', 'eta_R', 'eta_C'})
            % some parameters are defined as 1 - par, for these, the factor is
            adj_factor = (1-adj_factor)/val_to_adjust + adj_factor;
        end
    par_pf_changed.(par_to_adjust) = val_to_adjust * adj_factor; % adjust single parameter
    
    else
        % use absolute adjustment, leading to semi-elasticities
        adj_level = 0.01;
        
        if ismember(par_to_adjust, {'mu', 'eta_M', 'eta_R', 'eta_C'})
            % some parameters are defined as 1 - par
            par_pf_changed.(par_to_adjust) = val_to_adjust - ...
                adj_level;
        else
            par_pf_changed.(par_to_adjust) = val_to_adjust + adj_level;
        
        end
    end
    
    ProdFun = KrusellRobotsEach(par_pf_changed); % initialize production function object

    % Economy object
    econ.pf = ProdFun; % update economy object with adjusted production function object
    
    moments_change = compute_moments(econ, inp, lb, ub, q_B_0, q_B_1, droprange,w0_range,wf_range,prefix,suffix,dir_out,par_cal, par_pf_changed);
    moments = (table2array(moments_change) - table2array(moments_baseline))./table2array(moments_baseline);
    moments = array2table(moments);
    moments.Properties.VariableNames = moments_baseline.Properties.VariableNames;
end

function [results, results_percentiles] = simulate(econ,inp,lb,ub, q_B_0,droprange, w0_range,wf_range,prefix, suffix, dir_out, par, par_pf)
    sol = inp;
    stats = struct();
    percentiles = struct();
    results = table();
    results_percentiles = table();

    for i = 1:length(droprange)
        % change price of robots and recompute
        drop = droprange(i);
        econ.q_B = q_B_0 * (1-drop);
        
        fprintf('\n=================================================\n')
        fprintf('%s, price is %f\n', suffix, econ.q_B)
        fprintf('=================================================\n')
        
        if econ.fix_other_capital
            lb.K_E = inp.K_E;
            lb.K_S = inp.K_S;
            
            ub.K_E = inp.K_E;
            ub.K_S = inp.K_S;
        end
        
        tic
        [sol,fval,exitflag,output] = econ.compute_equilibrium(sol,lb, ...
                                                          ub);
        toc
        
        [stats, percentiles] = compute_stats(econ, sol, par, par_pf);
        stats.pricedrop = drop;
        stats.exitflag = exitflag;
        stats.fval = fval;
        results = [results; struct2table(stats)];
        results_percentiles = [results_percentiles; struct2table(percentiles)];
    end

    results_struct = table2struct(results);
    results_percentiles_struct = table2struct(results_percentiles);
    tools.struct2csv(results_struct,[dir_out,prefix,'_',suffix,'.csv'])
    tools.struct2csv(results_percentiles_struct,[dir_out,prefix,'_',suffix,'_percentiles.csv'])
end

function [stats, percentiles_strct] = compute_stats(econ, sol, par, par_pf)
    stats = sol;
    % assign prod fun pars
    stats.alpha = par_pf.alpha;
    stats.rho = par_pf.rho;
    stats.sigma = par_pf.sigma;
    stats.nu = par_pf.scale_L;
    stats.kappa_M = par_pf.kappa_M;
    stats.kappa_R = par_pf.kappa_R;
    stats.kappa_C = par_pf.kappa_C;
    stats.kappa_MR = par_pf.nu;
    stats.gamma_E = par_pf.lambda;
    stats.gamma_C = 1 - par_pf.mu;
    stats.gamma_B_M = 1 - par_pf.eta_M;
    stats.gamma_B_R = 1 - par_pf.eta_R;
    stats.gamma_B_C = 1 - par_pf.eta_C;
    stats.gamma_M = par_pf.upsilon;
    stats.A = par_pf.A;
    
    % assign other results
    stats.K_B = sol.K_B_M + sol.K_B_R + sol.K_B_C;
    stats.q_B = econ.q_B;
    stats.q_E = econ.q_E;
    stats.q_S = econ.q_S;

    [stats.Y_M,stats.Y_R,stats.Y_C,stats.Y_B_M,stats.Y_B_R,stats.Y_B_C,stats.Y_E,stats.Y_S]=econ.factor_prices(sol,par);
    stats.Y = econ.pf.Y(sol,par);

    w_lb = 1e-6;
    w_ub = 1e4;

    stats.part = econ.dist.mass_participants(stats.Y_M,stats.Y_R,stats.Y_C,w_lb,w_ub,par);
    stats.emp_M = econ.dist.emp_M(stats.Y_M,stats.Y_R,stats.Y_C,w_lb,w_ub,stats.part,par);
    stats.emp_R = econ.dist.emp_R(stats.Y_M,stats.Y_R,stats.Y_C,w_lb,w_ub,stats.part,par);
    stats.emp_C = econ.dist.emp_C(stats.Y_M,stats.Y_R,stats.Y_C,w_lb,w_ub,stats.part,par);

    stats.avg_wage_M = econ.dist.average_wage_M(stats.Y_M,stats.Y_R,stats.Y_C,w_lb,w_ub,stats.part,par);
    stats.avg_wage_R = econ.dist.average_wage_R(stats.Y_M,stats.Y_R,stats.Y_C,w_lb,w_ub,stats.part,par);
    stats.avg_wage_C = econ.dist.average_wage_C(stats.Y_M,stats.Y_R,stats.Y_C,w_lb,w_ub,stats.part,par);

    % welfare
    derived_vars = econ.compute_vars_from_input(sol);
    welfare = econ.welfare(derived_vars);
    stats.welfare_part_M = econ.welfare_participants_M(derived_vars);
    stats.welfare_part_R = econ.welfare_participants_R(derived_vars);
    stats.welfare_part_C = econ.welfare_participants_C(derived_vars);

    stats.welfare_non_part = econ.welfare_non_participants(derived_vars);
    stats.welfare_part = econ.welfare_participants(derived_vars);
    stats.welfare = econ.welfare(derived_vars);

    % tax revenues
    stats.tax_rev_inc = derived_vars.tax_rev_inc;
    stats.exp_transfer = derived_vars.exp_transfer;
    stats.tax_rev_B = derived_vars.tax_rev_B;
    stats.tax_rev_E = derived_vars.tax_rev_E;
    stats.tax_rev_S = derived_vars.tax_rev_S;
    stats.tax_rev_total = derived_vars.tax_rev_total;
    stats.excess_revenue = derived_vars.excess_revenue;

    % percentiles
    handle_h = @(w) 1/stats.part .* econ.dist.h(w,stats.Y_M, ...
                                                stats.Y_R,stats.Y_C,par);
    pct = 0.1:0.1:0.9;
    [percentiles, nodes, weights] = ...
        tools.Integration.compute_percentiles(handle_h,pct,1,200,100);

    % participation rates
    stats.part_M = econ.part_M(sol,par);
    stats.part_R = econ.part_R(sol,par);
    stats.part_C = econ.part_C(sol,par);

    % population shares
    integrand = @(w) econ.dist.dens_skill.f_M(w,stats.Y_M,stats.Y_R,stats.Y_C,par);
    stats.pop_M = tools.Integration.integral_legendre(integrand,w_lb,w_ub,econ.num_nodes_f,econ.nodes_f,econ.weights_f);

    integrand = @(w) econ.dist.dens_skill.f_R(w,stats.Y_M,stats.Y_R,stats.Y_C,par);
    stats.pop_R = tools.Integration.integral_legendre(integrand,w_lb,w_ub,econ.num_nodes_f,econ.nodes_f,econ.weights_f);

    integrand = @(w) econ.dist.dens_skill.f_C(w,stats.Y_M,stats.Y_R,stats.Y_C,par);
    stats.pop_C = tools.Integration.integral_legendre(integrand,w_lb,w_ub,econ.num_nodes_f,econ.nodes_f,econ.weights_f);

    percentiles_strct.percentile = pct';
    percentiles_strct.value = percentiles;
    percentiles_strct.q_B = econ.q_B * ones(size(pct'));
    percentiles_strct.q_E = econ.q_E * ones(size(pct'));
    percentiles_strct.q_S = econ.q_S * ones(size(pct'));
    
    % income shares
    stats.labor_income_share = econ.labor_income_share(sol,par);
    stats.labor_income = econ.labor_income(sol,par);
    stats.M_share_in_labor_inc = (sol.L_M * stats.Y_M)/stats.labor_income;
    stats.R_share_in_labor_inc = (sol.L_R * stats.Y_R)/stats.labor_income;
    stats.robot_share_capital = (stats.q_B * stats.K_B)/(stats.q_B ...
                                                      * stats.K_B + ...
                                                      stats.q_E * stats.K_E + stats.q_S * stats.K_S);

end

function moments = compute_moments(econ, inp, lb, ub, q_B_0,q_B_1,drop_range,w0_range,wf_range,prefix,suffix,dir_out,par_cal, par_pf)
    [results, results_percentiles] = simulate(econ, inp, lb, ub, q_B_0,drop_range,w0_range,wf_range,prefix,suffix,dir_out,par_cal, par_pf);
    percentiles_0 = results_percentiles(results_percentiles.q_B==q_B_0,:).value;
    percentiles_1 = results_percentiles(abs(results_percentiles.q_B-q_B_1)<1e-6,:).value;

    results_0 = results(results.q_B==q_B_0,:);
    results_1 = results(abs(results.q_B-q_B_1)<1e-6,:);
    
    percentiles_change = percentiles_1 - ...
        percentiles_0;

    K_B_0 = results_0.K_B;
    K_B_1 = results_1.K_B;

    K_B_rel_change = (K_B_1- K_B_0)/K_B_0;
    q_B_rel_change = (q_B_1 - q_B_0)/q_B_0;
    
    
    dir.generated = '../../../../data/generated/r/';
    dec_eps = tools.csv2struct([dir.generated,'dec_eps_AR2018.csv']);
    dec_eps = dec_eps.eps;

    targets.eps_w_K.val = dec_eps;
    
    dec_eps_results = percentiles_change./ ...
         percentiles_0./K_B_rel_change;
    
    bar([targets.eps_w_K.val,dec_eps_results])
    legend('data','model')
    
    % construct table with combined moments
    moments = results(1,:);
    moments.dec_eps_1 = dec_eps_results(1);
    moments.dec_eps_2 = dec_eps_results(2);
    moments.dec_eps_3 = dec_eps_results(3);
    moments.dec_eps_4 = dec_eps_results(4);
    moments.dec_eps_5 = dec_eps_results(5);
    moments.dec_eps_6 = dec_eps_results(6);
    moments.dec_eps_7 = dec_eps_results(7);
    moments.dec_eps_8 = dec_eps_results(8);
    moments.dec_eps_9 = dec_eps_results(9);
    moments.price_elast_robots = K_B_rel_change/q_B_rel_change;
    moments.part_change_M = (results_1.part_M - results_0.part_M)*100;
    moments.part_change_R = (results_1.part_R - results_0.part_R)*100;
    moments.part_change_C = (results_1.part_C - results_0.part_C)*100;
end